HIP-2008-03/TH 
1 February, 2008 
revised 1 July, 2008 



AN IMPROVED GLOBAL ANALYSIS OF NUCLEAR PARTON 
DISTRIBUTION FUNCTIONS INCLUDING RHIC DATA 

Kari J. Eskola'^'^ll], Hannu Paukkunen'^''^|^ and Carlos A. Salgado*^!! 



^Department of Physics, P.O. Box 35, FI-4OOI4 University of Jyvdskyld, Finland 

^Helsinki Institute of Physics, P.O. Box 64, FI-OOOI4 University of Helsinki, Finland 

Departamento de Fisica de Particulas and IGFAE, Universidade de Santiago de 

Compostela, Spain 



Abstract 

We present an improved leading-order global DGLAP analysis of nuclear parton 
distribution functions (nPDFs), supplementing the traditionally used data from deep 
inelastic lepton-nucleus scattering and Drell-Yan dilepton production in proton-nucleus 
collisions, with inclusive high-p^ hadron production data measured at RHIC in d-|-Au 
collisions. With the help of an extended definition of the function, we now can 
more efficiently exploit the constraints the different data sets offer, for gluon shadow- 
ing in particular, and account for the overall data normalization uncertainties during 
the automated minimization. The very good simultaneous fit to the nuclear hard 
process data used demonstrates the feasibility of a universal set of nPDFs, but also 
limitations become visible. The high-pT forward-rapidity hadron data of BRAHMS 
add a new crucial constraint into the analysis by offering a direct probe for the nu- 
clear gluon distributions - a sector in the nPDFs which has traditionally been very 
badly constrained. We obtain a strikingly stronger gluon shadowing than what has 
been estimated in previous global analyses. The obtained nPDFs are released as a 
parametrization called EPS08. 
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1 Introduction 



With collider energies presently reached at BNL-RHIC, hard processes have become 
more and more important as diagnostic tools in the phenomenology of heavy ion colli- 
sions, QCD matter and QCD dynamics. The soon starting LHC heavy-ion program will 
emphasize the role of hard processes even further, by extending the kinematical range 
probed in the longitudinal momentum fraction x and in the process virtuality scale 

by several orders of magnitude with respect to the presently accessible ones. The 
existence of well-constrained up-to-date nuclear parton distribution functions (nPDFs) 
will thus be essential for the correct interpretation of the data. 

Sets of collinearly factorized universal nPDFs, which are obtained in global pertur- 
bative QCD analyses paralleling those for the free proton, are available [U, E], [3], HI El 
|6],[7], see also Ref. [8]. These nPDFs are typically obtained by first parametrizing the 
nuclear corrections for each parton flavour relative to a known set of the free proton 
PDFs, and imposing constraints from sum rules at a chosen initial scale Q^. A best 
fit to nuclear hard-process data from deep inelastic lepton-nucleus scattering (DIS) 
and the Drell-Yan (DY) process in proton-nucleus collisions is obtained by an iterative 
procedure which involves the DGLAP evolution [9] of the (absolute) nPDFs. The best 
global fit then fixes the initial nuclear corrections. 

Since the first of the nPDF sets, EKS98 [H [2] , the procedure has been improved 
by performing the analysis at next-to-leading order (NLO) [HI [7] and by making un- 
certainty estimates [3], [5l [7] in analogy with the free proton case. In spite of such 
important progress, however, new data sets of relevance, which would more directly 
constrain the nuclear gluons in particular, have not been included in these ten years. 
The purpose of this paper is to make progress precisely in this respect, by including 
the data from inclusive high-pT hadron production in d-l-Au collisions at RHIC in the 
global analysis of the nPDFs for the first time. 

The most serious difficulty in the global DGLAP analyses of nPDFs has tradition- 
ally been the lack of experimental data which would impose stringent enough con- 
straints for the nuclear gluon distributions. The extraction of the nPDFs would be 
cleanest in DIS but basically no high-precision data are at hand in the perturbative re- 
gion > 1 GeV^ at a; < 0.01 there. This deficiency translates into a bad determination 
of the nPDFs, gluons in particular, in a region, where the nuclear effects are sizeable 
and which will be frequently accessed at the LHC. For the analysis of hard processes 
taking place at midrapidities at RHIC the situation has been better, as the available 
DIS and DY data constrain the nPDFs in most of the kinematic range probed. In 
the forward-rapidity domain, however, the hard processes are sensitive to nPDFs at 
smaller values of x than what is presently constrained by DIS and DY: These RHIC 
data offer a possibility for independent further constraints of the nuclear gluon sector 
in the global analysis. 

As discussed in previous works, |3l [10], the suppression in high-p^ hadron pro- 
duction at forward rapidities in d-l-Au collisions relative to p+p collisions, measured 
by BRAHMS [TT], would seem to suggest stronger gluon shadowing than, e.g., in the 
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EKS98 set. This suppression has been also proposed as a signal of parton saturation 
being reached at RHIC [I2], so that the compatibility of the measured suppression with 
DIS and DY data sets within the DGLAP framework is a question of special relevance 
from the QCD parton dynamics point of view as well. 

In this paper, we shall demonstrate for the first time, that a global fit of a very 
good quality can indeed be obtained by simultaneously accommodating the DIS, DY 
and high-pT RHIC data in the leading-order (LO) DGLAP framework, i.e. that a 
relevant new set of universal, process- independent, coUinearly factorized nPDFs can 
indeed be extracted, and that the gluon shadowing obtained at the smallest values of 
X is indeed stronger than in previous global fits. Also limitations and uncertainties 
remaining in the analysis are discussed in light of the results obtained. An important 
new feature which we introduce in the global x^-analysis here - borrowing it from 
the free proton analyses [13] - is the treatment of data normalization errors given by 
the RHIC experiments. In particular, we demonstrate that only by accounting for 
these systematic errors, a meaningful comparison with the RHIC high-p^ pion data 
O [151 [16] can be done. 

Parametrizing the obtained nuclear effects for each parton flavour in x, and A 
and making a simple fast computer code available for public use has proven to be a 
working idea in the past. Analogously with our previous EKS98 set [B [2], we now 
release a new set of nPDFS called EPS08, which is available at |17j . 

The rest of this paper is organized as follows: In Sec. [2] we present the framework 
and analysis method, introducing the functional forms used and, in particular, the 
improvements in the fitting procedure. In Sec. [3] we present the results from the 
global fit, the new set of nPDFs, and the comparison with experimental data. The 
effects of the normalization-error treatment are demonstrated. In Sec. HJwe comment 
on the strong gluon shadowing solution found. Conclusions and outlook are presented 
in Sec. |5l 



2 Framework and analysis method 
2.1 Definition of nPDFs 

The DGLAP framework we use is essentially the same as in our previous global fits 
of nPDFs [H [21 [3]- For each parton flavour i, we define the nPDFs ff^{x, Q^) as the 
PDFs of protons bound to a nucleus of mass number A, 

/^(x,Q^)^i?f(x,g^)/f™Q^^^(x,Q^), (1) 

where f^"^^^^^^{x,Q'^) is obtained from the latest LO CTEQ set of the free proton 
PDFs [18], and Rf{x,Q^) is the nuclear modification factor for this parton flavour. 
We also assume that PDFs of bound neutrons can be obtained on the basis of isospin 
symmetry. For instance, the total u quark PDF in a nucleus A with Z protons then 
becomes 

ua{x, Q') = Zf^ix, Q') + {A- Z)f^{x, Q'). 
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The nuclear effects of Deuterium [A = 2) and tliose in tlie cumulative region x > 1 
are neglected. We do not discuss the dependence of the nPDFs on transverse location 
inside the nucleus (the impact parameter dependence) [101 [H] here, either, only the 
average nuclear effects are considered. 

We parametrize the nuclear modifications Rf at an initial scale Qq = 1.69 GeV^, 
which conveniently matches the lowest scale and the charm-quark mass threshold in 
the CTEQ6L1 set. The valence quark modifications are constrained by baryon number 
conservation and the gluon modifications by momentum conservation. At higher scales 

> Qq the nPDFs are then obtained by solving the conventional DGLAP equations 
[9] at LO numerically, applying the fast solution method introduced in [20] . 

In order to reduce the amount of fitting parameters, and obtain a converging well- 
constrained fit, we have to assume initially flavour-independent nuclear effects for va- 
lence quarks, sea quarks and gluons, i.e. we parametrize only three different functions: 
Ry{x, Ql) for all valence quarks, Rg^x, Ql) for all sea quarks, and Rq{x, Ql) for gluons. 
In the DGLAP evolution, each flavour is considered individually so that at > Qq 
the modifications may in principle depend on parton flavour. 

2.2 Fitting functions and parameters 

Choosing a suitable functional form for the input nuclear modification factors Rf{x, Qq) 
is among the most troublesome and crucial issues in the global nPDF analyses. On 
one hand the fit functions must be flexible enough, i. e. there should appear sufficiently 
many parameters so that all the relevant features suggested by the data can be caught. 
On the other hand, too large a number of parameters easily leads to badly converging 
fits. Thus the number of parameters always is a compromise between the flexibility of 
the fit and the feasibility of the x^-analysis. 

Up to now, gluon shadowing in the small-x region x < 10^^ has been constrained 
only by the momentum sum rule, and it has been essentially dictated by the assumed 
form of the fit function. Also the statistical error bars computed then reflect the uncer- 
tainties within the fit function chosen. Thus, as discussed in [3], there has been a large 
uncontrolled error in the gluon shadowing. The inclusion of the high-p^ hadron data 
from RHIC at forward rapidities, however, now provides important further constraints 
for the gluon shadowing region, which will be exploited in the present work. 

In order to take into account the RHIC high-p^ forward-rapidity hadron data, which 
suggest a stronger gluon shadowing than obtained in previous global analyses (see the 
discussion in [5]), we need to modify the parametrization of the shadowing region. 
In our past works, we assumed a saturation of the modifications Rf{x,Ql) such that 
Rq ~ Rg — > const at X ^ 0. This assumption is relaxed in the present analysis, and 
we introduce a power-law behaviour of the nuclear modification factors at small x as 
follows 

i?f(x,g2)-^%"-^, «^>o. (2) 

This can be motivated by the (approximate) power-law behaviour ~ x^^'^^" of the free 
proton PDFs at small-x, assuming that the nPDFs share this same gross feature but 



3 



This means that Rf^ 



~ X 



{Pfr 



l) 



as 



with a different power Pfrce > -Pbound- 
x-^0. 

With this assumption as a guide, we parametrize the initial nuclear modifications 
Ry, Rg and Rq in three pieces as illustrated in Fig. [1) Rf{x) at small values of x, 
below the antishadowing maximum, x < x'^] R^i^) from the antishadowing maximum 
to the EMC minimum, x^ < x < xf; and Rf{x) in the large-x Fermi- motion region, 
X > ; 



Rf{x) 
Riix) 

Rfi 



Cq + (cf + c^x° )[exp{—x/xf)—exp{—x^/xf)], x < x^ 



x) 



°0 



afx + a^x^ + a^x^. 



Xf> 



xf < X < x:5 



X„ < X < 1. 



(3) 



Some of the parameters above are eliminated by matching the different pieces smoothly 
together: we require continuity of the fit functions and zero first-derivatives at the an- 
tishadowing maximum x^ and at the EMC minimum xf . The required behavior ([2]) 
It is convenient to express the fit functions in terms of the following 8 param- 
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Figure 1: An illustration of the smoothly matched fit functions Rf{x) and the role of the 
parameters x^, y^, and A^. The superscripts A have been suppressed in the figure. 



In principle, all parameters above are different from one nucleus to another, hence 
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the superscript A. For the A dependence we assume a simple power law, 



zf = zt-'i^r% (4) 

^ref 

where Zi = Xs, Xa, Va ■ ■ ■, and choose the reference nucleus to be Carbon, Arcf = 12. 

With momentum and baryon number sum rules, we can fix for gluons and 
valence quarks, and thus reduce the total number of parameters to 44. This is still far 
too many for a convergent x^-minimization with the nuclear data constraints available. 
To proceed, additional assumptions need to be introduced - how the unconstrained 
regions of the phase space are handled, and how the number of final fit parameters 
is reduced down to 15 by fixing those parameters which cannot be constrained, is 
explained in Sec. 13.11 

2.3 Data sets 

The experimental data, providing the nonperturbative input for the nPDFs, in our 
present analysis covers three types of hard processes involving nuclei: In addition to 
the data from lepton-nucleus deep inelastic scattering and proton-nucleus Drell-Yan 
dilepton production, as a new ingredient we include the data from inclusive high-pr 
hadron production in minimum-bias d+Au collisions from the BRAHMS, PHENIX 
and STAR collaborations at RHIC. In total, we have over 600 data points covering 13 
nuclei from Helium up to Lead. Table [1] summarizes the data sets in our analysis. 

The (minimum bias) DIS and DY data that we utilize are available as ratios of 
differential cross sections between a nucleus A and a reference nucleus. We denote the 
cross section ratios computed against deuterium as 

^ = Jd^^JdQ^ - ^^^^^'^ ^' ^ = Idal^/dM^dx ' 

(5) 

where x refers to the momentum fraction and to the photon virtuality for DIS, and 
to the invariant mass of the lepton pair and x to either xi or X2 for DY. The scales 

and also define our factorization scales, and we consider only those data points 
which lie above our initial scale: M^,Q^ > 1.69 GeV^. 

The inclusive hadron production data at RHIC comes as the nuclear modification 
factor -RdAu; the ratio between the invariant yields in d-l-Au and p+p collisions, 

^ 1 d'^m^^/dprdr] „,i,^bias ^d^a^'-ydprdri 

(A^coii) d^NPP/dpTdT] d^aPP/dprdr] ' ^ ' 

where pt and r] denote the transverse momentum and the pseudorapidity of the ob- 
served hadron, and (A'^coii) is the estimated average number of inelastic binary nucleon- 
nucleon collisions in a centrality class studied. The last equality holds for the minimum 
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Table 1: The data used in this analysis. The mass numbers are indicated in parentheses 

and the number of data points refers to those falling within our kinematical cuts, Q^, > 

1.69 GeV^ for DIS and DY, and px > 2GeV for hadron production at RHIC. The quoted 

values correspond to the unweighted contributions of each data set. 
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bias case which we are interested in here. According to the QCD factorization theorem, 
the inclusive hadron production cross sections can be computed as 

ijkl 

where ff^^Q"^) are the input nPDFs, cr*-'~**^'(Q^) is the pQCD matrix element squared, 
and Dk-th+xiQ}) denotes the fragmentation functions for which we use the KKP 
parametrization [21]. Our choice for the fragmentation scale is the hadronic trans- 
verse momentum, Q'j = and for the factorization and renormalization scale we 
take the corresponding partonic transverse momentum. Notice here that all the scale 
choices above are simplifications in the sense that they could be left as additional fit pa- 
rameters, and that the KKP fragmentation functions do not distinguish negatively and 
positively charged hadrons. Such details, however, are beyond the scope of the present 
analysis. A more detailed discussion of how the needed differential cross sections ([7]) 
are calculated in practice in LO can be found in |22j . 

We should emphasize that the following choices are made with the inclusion of the 
RHIC data: 

• Since we do not discuss the impact parameter dependent nuclear effects here, we 
consider only minimum bias data for -RdAu- 

• From PHENIX and STAR, we systematically include only the pion production 
data. This is because the 'Cronin-type' enhancement seems to be much larger for 
baryons than for mesons [11] , signalling of the fact that still at V^nn ~ ^00 GeV 
there might be a significant component of nonperturbative baryon-number trans- 
port from the beam particles. We do, however, include the BRAHMS data for 
negatively charged hadrons with the assumption that the antiproton content in 
this data sample is negligible. 

• Perturbative QCD calculations of inclusive hadron production, which are per- 
formed strictly in the framework of coUinear factorization (without any intrinsic 
transverse momentum), show similar features both in LO [22] and in NLO [231 [23]: 
Below pj- ~ a few GeV the computed cross sections for p + collisions start 
to overshoot the data. Hence, since we should not push the nuclear case too far 
either, we choose to include only the region > 2 GeV of the RHIC data on 
inclusive hadron production in this analysis. More discussion on this choice will 
follow later. 

2.4 Modified 

The established way of fitting a set of parameters {z} of the PDFs against a large 
number of experimental data, is the minimization of the global function. In its 
simplest form the global x^, the goodness parameter of the fit obtained, is defined by 

X\{z]) ^Y.^l{{z]\ (8) 

N 
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where labels the experimental data sets, and 



D,-Ti{{z}) 



(9) 



where Di, ai, and Ti{{z}) denote the value of a single data point, its measurement 
uncertainty, and the corresponding theoretical value which depends on the parameters 
{z} of PDFs. 

For the nPDF analyses involving only DIS and DY data, the simple form of 
above has been sufficient, but for our current purposes a more general definition for 
the x^, introduced in [13], is needed: 



where wn is a weight factor chosen separately for each data set, 0"^°™ is the relative 
uncertainty in the overall normalization reported by the experiment, and is the 
optimized value of the overall normalization for the data set, corresponding to each 
parameter set {z}. The reasons for the necessity of such redefinition are the following: 

1. By making the weight factor wn larger than 1, we can emphasize by hand the 
importance of those data sets which contain definite physics content - such as 
constraints for small-x gluons - but whose number of data points is small. With 
a default value = 1, such data sets would have a negligible contribution to 
the overall the valuable constraints they offer would escape unnoticed. 

2. In addition to the point-to-point statistical and systematic errors, certain data 
sets have a significant common normalization uncertainty a]^''™ for all data points 
within the set. Even if this normalization uncertainty is large, the shape of 
the distribution formed by the data points may be a valuable constraint for 
the nPDFs. We introduce for each data set a normalization factor /^r G [1 — 
^norm^ 1 + 0"]^°™] which multiplies all the experimental values within the set N. In 
connection with /at, there is an additional "penalty" factor ( ^no^m )^ which is the 
larger the more deviates from unity — this accounts for the fact that having 
/tv = 1 is anyway the experiment's best estimate for normalization. The actual 
value for /tv is determined from the requirement that XnH^}) for each data set 
is at minimum. 

The motivation for both modifications in the definition discussed above comes 
mainly from adding the RHIC data for the nuclear modification factor -RdAu of Eq. 
into the analysis. First, the BRAHMS data set for forward direction (77 ~ 2. ..3, espe- 
cially with our choice > 2 GeV) has only a very few data points. These would not 




(10) 



TV 




(11) 
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have much effect in the global without being artificially emphasized. Second, the 
average number of inelastic binary nucleon-nucleon collisions (A'^cou) in d + Au collision 
is derived from a simulation of the experiment with the Glauber model as an input, 
which gives rise to a significant model-dependent normalization uncertainty in -RdAu- 

The amount of DIS data overwhelms that of the DY data. To improve upon this 
balance, we weight the FNAL-E772 DY data set hy wn = 10. This improves the deter- 
mination of the relative importance between the valence and sea quarks at intermediate 
values of x. For the NMC data set for Rp^ we give a weight wn = 5 in order to better 
ensure a good fit for the Carbon nucleus, which is used as reference in the analysis. 
The NMC 96 data on the Q^-dependence of F^^'/F^ is weighted by = 10 but only 
for the three lowest values of x - the three upper panels in Fig. [8] below - to help 
constraining the gluon distribution at x < 0.02 via the DGLAP evolution. Finally, the 
few points of the BRAHMS data that we include in our analysis, are weighted by a 
large factor w;7v = 40 in order to account for the constraints this data set gives for the 
gluon distribution. All these weights are summarized in Table [1] above. 

3 Results 

3.1 Final parameters 

In order to reach a well converging (well constrained) global fit, where none of the fit 
parameters are drifting to their limits, we are forced to reduce the total number of free 
parameters down to the following 15: 

• Valence quark modification 

The DIS data constrain the modification Qq) in the a; > 0.1 region rather 

well, and altogether 8 parameters Xa, Va, Vya^ ^e, ^e, PAe, ^2, Pb2 were left free. 

• Sea quark modification 

The DIS and DY data probe the sea quarks in the region 0.01 < x < 0.1, and 
5 parameters, a, pa, Xa, Va, Vya^ controlling this region in Rg^x.QD, were left 
free. The region x^ 0.3 is, however, not constrained by any present experimental 
data. We assume for simplicity a smooth behavior in this region, without an 
EMC effect. 

• Gluon modification 

The gluon modification Rq{x, Qq) at small x is now directly constrained by the in- 
clusive hadron production data from RHIC. Indirectly the gluons are constrained 
by the Q^-evolution effects in the sea quark sector, reflected by the DIS and DY 
data. In spite of the new constraints, we were still able to leave only 2 parameters, 
Ua and controlling the antishadowing peak height, free. We assume a similar 
EMC-effect for gluons as there is for valence quarks, guided by the shape of the 
preliminary PHENIX data for inclusive photon production in Au+Au collisions 
at ~ 6 GeV [33] and also by the PHENIX data for inclusive pion production 
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at ?7 = and ~ 6 GeV [15] (see Fig. [12] ahead - we have checked that indeed 
some sensitivity to gluon PDFs persists even at these values of pr)- 

The global fit with these 15 parameters was then performed by minimizing the 
function defined in Eqs. ^ and ^ with the MINUIT [M] routine from the CERN 
Program Library. Table [2] summarizes the parameter values obtained as well as the 
fixed parameters. The goodness of the fit is characterized by /N = 0.71, where 
is computed with no extra weights, wn = 1, but with the optimized normalization 
factors /tv included, and = 627 is the total number of data points. The contribution 
from each data set to this can be read off from Table [1] The corresponding nuclear 
modifications for selected nuclei at our initial scale Ql = 1.69 GeV^ are shown in Fig. El 
Table [3] shows the contribution from different types of hard processes to the unweighted 

and the comparison with our previous analysis in Ref. [3]. The /N obtained in 
earlier global analyses can be found in Table 3 of Ref. [3] . 
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Table 2: List of all parameters defining the modifications Ry, Rg and Rq through Eq. ([3]) 
at our initial scale Qq = 1.69 GeV^. The parameters a, x^, Xq, Xg, Va-, Ae, 62 and (3 are for 
the reference nucleus A = 12, and the powers pi define their ^-dependence as in Eq. For 
valence quarks and gluons, the baryon number and momentum sum rules fix the parameters 
a"^ for each nucleus A separately, in which case the powers pa are not used. The location 
and height of the EMC minimum of Rq was fixed to that of Ry. The parameters left free 
for the minimization procedure are shown in bold face. 
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Data type 


Data points 


Xepsos 


Xekps 


Deep Inelastic 


484 


344.2 


337.4 


Drell-Yan 


92 


77.1 


84.3 


Hadron production 


51 


26.9 


28.0 


Total 


627 


448.3 


449.6 



Table 3: Contributions of various data types to the total unweighted in our previous 
work [3] (EKPS) and in this work (EPS08). 



3.2 Comparison with data 

The comparison of our results with the experimental data used in the global fit is 
presented in Figs. [3] - E] for DIS, in Figs. [9] - [ID] for DY, and in Figs. HH - [12 for 
the RHIC data. In all figures, the open (red) symbols denote the experimental data. 
The error bars for the DIS and DY cases correspond to the point-to-point statistical 
and systematic errors added in quadrature, while for the RHIC data the statistical and 
systematic errors are shown separately. 

In comparison with our previous work, Ref. [5], the Copper and Lithium data in 
Figs. IHEl have now been added into the analysis. In Fig. [31 we see that shadowing for 
heavy nuclei (Sn and Pb) has now gotten stronger, and in Fig.[7|that we now reproduce 
the largest-x DIS data better than before. The smallest-x panel of Fig. [8] is one of the 
key issues in this paper and the obtained slopes will be separately commented in 
Sec. [H below. Figure [HI shows that for the DY ratios the A systematics at the smallest 
values of X2 have been improved: now also the Tungsten data are reproduced well. 
This improvement is reflected also in the large-xi part of the W/Be ratio in Fig. [TOl 

The data set that plays a major role in constraining the gluon modifications in 
the present analysis, is the inclusive negatively-charged hadron production at forward 
direction (77 = 2.2 and 77 = 3.2) measured by the BRAHMS collaboration at RHIC, 
shown in Fig. [TH For the data sample we include in the global fit, > 2 GeV, the 
optimized normalization factor is close to one, /at = 1.02. 

Figure [12] presents the comparison with the PHENIX and STAR measurements of 
inclusive pion production at midrapidity (77 ~ 0). The need of a treatment which ac- 
counts for normalization uncertainties is clearly demonstrated by this figure. Although 
all data sets agree within the given large uncertainties, the general trend in the STAR 
data is somewhat different from the PHENIX data, as can be seen in the uncorrected 
case (/at = 1), shown on the left-hand side of the figure. Taking into account the 
normalization uncertainties as provided by the modified definition of in Eqs. ( 1111) . 
a good fit with both PHENIX and STAR data sets becomes indeed possible - see the 
right-hand side of Fig. [T21 where the optimized normalization factors for the PHENIX 
data are = 1.04 and 1.07 and for the STAR data = 0.90. 
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X X 



Figure 2: The nuclear modification factors Ry, and Rq for C, Ca, Sn, and Pb at 
Ql = 1.69 GeV^ The DIS ratio R^^ is shown for comparison. 

4 Discussion 

Figure [13] shows a comparison of the nuclear effects in the average valence quark, 
average sea quark and gluon distributions at our initial scale Qq = 1.69 GeV^, as ob- 
tained for a Lead nucleus in the LO DGLAP analyses here (EPS08), in HKN07 [7], 
in nDS [B], and in our previous works EKPS [5] and EKS98 [2]. The figure demon- 
strates the fact that while the average effects in the valence quarks and in the mid-x 
region (0.01 ;5 x 0.2) of the sea quarks are relatively well under control, quite large 
uncertainties remain in the large-x {x^ 0.2) sea quark and gluon distributions. In par- 
ticular, the gluon shadowing which we obtain here on the basis of the BRAHMS data, 
is clearly stronger than in the previous global analyses. A strong gluon shadowing 
has been suggested before at least in the Glauber-Gribov framework [351 ES] (see also 
the review [27]) and also in the context of DGLAP evolution [3H], but not in a global 
DGLAP analysis where constraints from DIS, DY and RHIC hadron production data 
are simultaneously imposed. 

The obtained gluon shadowing reflects a compromise between a weaker gluon shad- 
owing suggested by the NMC 96 DIS data [27] (the first panels in Fig. [S]), and a 
stronger effect demanded by the BRAHMS data. Regarding these constraints, we note 
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Figure 3: The computed ratio Rp^{x,Q'^) vs. Rp^{x,Q'^) compared witli tlie NMC data 
|28| . Tlie open symbols are tlie data points with statistical and systematic errors added in 
quadrature, the filled ones are the corresponding results from this analysis. 



13 



0.85 



c 




\ 


O - 

• - 


A=6 

$ * ■ ■ 

-I * 

[] 


[ 

1 

1 

1. 


■ - 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 II 


) 




A=40 ,H ^ 

• T ? T 


- [] 


^ ^ o NMC95re 
□ NMC95 


'o 

1 

!! 


II 

[] - 
o - 




I 

c 


C) 

- 



10"' 10" lo' 10" 10" 10 



Figure 4: The calculated ratio Rp^{x,Q'^) compared with the NMC 95 (squares) [27] and 
the reanalysed NMC 95 (circles) data |26| . 
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NMC data (open circles). 



14 



\j.y r 



0.8 - 



0.7 
10 



- □ EMC 



■3 



10- 



10" 



Figure 6: The calculated ratio Rp^{x, Q"^) (filled squares) compared with the EMC [HI] data 
(open squares). 
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Figure 7: The calculated ratios Rp^{x,Q'^) (filled squares) for several nuclei compared with 
the SLAC data [25] (open squares). 
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Figure 8: The calculated scale evolution (solid black lines) of the ratio F2^'^/F^ compared 
with the NMC data [32] for several fixed values of x. 
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Figure 9: The computed R^Yi^2, M"^) (filled squares) as a function of X2 compared with the 
E772 data [2^ (open squares). 
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Figure 10: The computed R^yi^i^ ^'^) (filled squares) as a function of xi compared with 
the E866 data [30] (open squares) at four different bins of invariant mass M^. 
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Figure 11: The computed nuclear modification ratio -RdAu at forward rapidities (filled 
squares) for negatively-charged hadron production, compared with the BRAHMS data [11] 
(open squares). The error bars are the statistical uncertainties, and the shaded bands indi- 
cate the point-to-point systematic errors. The additional overall normalization uncertainty, 
is 5%, i.e. a^^"^ = 0.05 in Eq. (Ilip . The upper panels show the comparison without the 
normalization factor fj^. In the lower panels, the data have been multiplied by the optimized 
value /at = 1.02. 
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Figure 12: The computed i?dAu (filled symbols) at midrapidity (r/ = 0) for inclusive pion pro- 
duction compared with the PHENIX [141 115j and STAR [16] data (open symbols). The error 
bars are the statistical uncertainties, and the shaded bands indicate the point-to-point sys- 
tematic errors. The additional overall normalization uncertainties are 10% for the PHENIX 
data and 17% for the STAR data. The left panels show the comparison without the normal- 
ization factor /jv- In the right panels, from top to bottom, the data have been multiplied by 
the optimized values /at = 1.04, = 1-07 and = 0.90. 
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Figure 13: Comparison of the average valence and sea quark, and gluon modifications at 
Q2 = 1.69 GeV2 for Pb nucleus from LO global DGLAP analyses EKS98 [2], HKN07 [7], 
nDS [6], EKPS [3] and this work EPS08. 




X Q pj' 



Figure 14: The correlation between the gluon shadowing (left), the slopes of /F2 
(middle) and the forward- 77 -RdAu (right) when weights ■wat = 0, 40 and 150 are assigned 
for the BRAHMS data. In the right panel, the comparison is shown without an overall 
normalization factor, i.e. with /at = 1. 
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the following: 

First, the DIS data [27] show that the dependence of Rp^ at x ~ 0.01 in the 
region ~ 1 GeV^ is very weak. Given this, the small-x approximation of the DGLAP 
equations [391 HQ], 

91ogg2 - 2771 F^{x,Q') > )\ , U^J 

then indicates that gluon shadowing is restricted to be similar to what has been mea- 
sured for F2. 

Second, using Eq. (fT2l) above, and the fact that the logQ^-slope of / F2 mea- 
sured by NMC has been observed to be positive at x = 0.125 - see the first panel in 
Fig. [8] - we deduce that 



i?g(2x,Q2) 



sO.0125 -"-^2 



R%{x,Q^) 



(13) 

a;Ri0.0125. 



This indicates that the A dependence of gluon shadowing is weaker than that of F2. 

Third, the A dependence of gluon shadowing must still be strong enough in order 
to reproduce the BRAHMS data. The flattening of the Q^-slope seen in the first panel 
of Fig. [8] indicates that the gluon shadowing now obtained - the A dependence of the 
gluon modifications in particular - is already so strong that it is in the brink of violating 
the condition f|T3|) . In this sense the gluon shadowing in the present global fit is the 
strongest possible one which is still in agreement with the DIS data. 

Fourth, a balance between the constraints that the NMC 96 and BRAHMS data 
offer for the gluon shadowing, is obtained by assigning suitable relative weights, see 
Table [TJ Since these data sets drive the fit to opposite directions, the resulting gluon 
shadowing obviously depends on the weights introduced. To demonstrate this sensitiv- 
ity, we have repeated the analysis by varying the BRAHMS data weights as follows: By 
setting wjv = 0, we remove this data set from the analysis. Alternatively, by assigning 
a very large weight, wn = 150, to the BRAHMS data, we clearly overemphasize its 
importance. In both cases, the smallest-x NMC 96 data set weight is kept unchanged 
{wn = 10, see Table [1]). The overall fits obtained in these extreme cases remain very 
good, giving x^/N = 0.72 and 0.73, correspondingly. Figure [H] (left panel) shows the 
resulting gluon modifications in each case, along with a comparison to the smallest-x 
NMC 96 data (middle panel) and the BRAHMS data (right panel). The figure clearly 
demonstrates how adding more weight to the BRAHMS data will eventually flip the 
sign of the computed slopes of F^^^/F^ — a phenomenon which on the basis of the 
systematics seen in the NMC 96 data would be an unwanted feature. With a weight 
factor Wn = 40 for the BRAHMS data (making the effective number of the BRAHMS 
data points the same as in the three smallest-x panels of the NMC 96 data), we reach 
the strongest possible gluon shadowing, and thus a fair agreement with the measured 
forward-rapidity -RdAu, without such a sign flip. 

Naively, in the RHIC hadron data, we could well expect that hadrons at fixed rj and 
Pt would dominantly come from partons of higher transverse momenta, qx ~ 1.5. ..2pT 
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and the same rapidity t] [22], whose production would mainly probe the nPDFs at 
momentum fractions X2 = + e"^^) ~ ~ 10~^, assuming 2 — 2 parton 

production kinematics, and taking px ^ 2 GeV and 1/2 ~ ''7 = 3.2. We notice, however, 
that the ratio -RdAu at = 2 GeV in Fig. [11] is considerably larger than the gluon 
shadowing we obtain at these values of x, see Fig. [13] (solid line). This is due to 
two reasons: First, as shown in Fig. 13 of [3], the DGLAP evolution from Qq to the 
few-GeV region increases the ratio Rq substantially. Second, the integration over the 
partonic qt and over the unobserved parton rapidity ?/2 causes a significant smearing 
of the x-range probed, especially towards larger x (see also Ref. [H] and Table 1 in 
|10j). Thus, the ratio -RdAu at forward 77 is in fact sensitive not only to nuclear gluon 
shadowing but also to gluon antishadowing - and antishadowing in turn amplifies when 
shadowing gets stronger. Therefore, even a large change in gluon shadowing induces 
only a moderate change in the computed ratio -RdAu? and a significant gluon shadowing 
is required in order to reproduce the BRAHMS data at > 2 GeV. 

As seen in Figs. [11] and [121 we have a good fit of the nuclear modification factor 
-RdAu at pt > 2 GeV. We cannot, however, reduce -RdAu by strengthening the gluon 
shadowing as much as the BRAHMS data below 2 GeV would require without violating 
the DIS data constraints. This is also one of the reasons for excluding the region < 2 
GeV of the RHIC data from this analysis. As explained above, we also are hesitant 
to push the nuclear case too far into the small-pT region, for we cannot reproduce 
the shape of the absolute pt spectra in p+p collisions well enough there, and for we 
do not consider impact-parameter dependence of nPDFs or a more detailed centrality 
selection here. 

As illustrated by Fig. [131 a saturation of shadowing at x ^ was assumed in 
previous global analyses. This assumption is now relaxed with the aim to study the 
strongest gluon shadowing allowed by present experimental data. It is worth empha- 
sizing that the behavior of the nuclear corrections at the smallest values of x (x < 10~^) 
is still largely determined by the assumed shape of the fit functions in Eq. 1^. This 
limitation is common to any global fit (for nPDFs as well as for the free proton PDFs) 
in those regions of phase space which are poorly or not at all constrained by the data. 

5 Conclusions 

We have improved the global analysis of nPDFs in two important ways: First, by taking 
the RHIC data into account in such analysis for the first time, we have extended the 
constrained x region down to 10^^. Second, we have improved the minimization 
procedure by introducing weighting of different data sets and by explicitly accounting 
for the overall normalization errors quoted by the experiments. 

One of the main goals of this paper is to study to what extent a strong gluon 
shadowing suggested by the BRAHMS data can be accommodated together with the 
DIS and DY data in a global analysis. We conclude that a simultaneous fit of DIS, DY 
and high-pT- {px > 2 GeV) hadron production data from RHIC at forward rapidities 
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(BRAHMS negative hadrons), is indeed possible within the DGLAP framework without 
invoking any new suppression mechanism. Thanks to the improved treatment of data 
normahzation errors, we obtain also a good agreement with the RHIC pion data at 
mid-rapidity (STAR and PHENIX). The very good quality of the global fit obtained 
suggests that a well-working universal set of nPDFs can be extracted in the framework 
of coUinear factorization. Within an improved global analysis, and with an emphasis 
on the RHIC forward-rapidity data, we obtain a stronger gluon shadowing than in the 
previous global nPDF analyses, see Fig. [131 These are the main new results of this 
paper. The LO nPDF set we have obtained (EPS08), i.e. a parametrization of the x, 
and A-dependent nuclear modifications relative to CTEQ6L1, is available at [17] 
for practical use. 

As discussed above, the amount of gluon shadowing obtained depends on the weight 
assigned to the BRAHMS data, and equally good overall fits can be obtained also when 
the weights are smaller and the resulting gluon shadowing is weaker. Until more data 
become available to resolve this problem, the nuclear gluon distributions suffer from 
considerable uncertainties. To estimate the effects of these uncertainties in the hard 
process cross sections one computes in LO, we recommend to use the current results, 
EPS08, in parallel with the previous LO results EKS98 [2], nDS [6] and HKN07 [7j. 

Regarding inclusive hadron production in nuclear and hadronic collisions, computed 
here in the coUinear factorization framework, we would like to emphasize that we have 
limited the study to the region pr > 2 GeV: Within the present global analysis, we 
cannot reproduce the sudden drop of the ratio -RdAu measured for negative hadrons 
by BRAHMS at px < 2 GeV at forward rapidities, see Fig. [HI or the very strong 
suppression of -RdAu measured by STAR for 7r° at r/ = 4 and pr < 2 GeV ^2] - a yet 
stronger gluon shadowing needed for this would clearly lead into a contradiction with 
the logQ^ slopes of F2^/F2 measured at DIS. More detailed work on fragmentation 
functions, impact parameter dependence of the nPDFs as well as further developments 
in the fit functions is required in order to make firmer conclusions on the applicability 
of the DGLAP-evolved universal nPDFs in this region. Regarding the fragmentation 
functions, we anticipate that considering a more detailed charged separation (see e.g. 
[l3l HU Us]) in hadron production would tend to increase the computed -RdAu rather 
than decrease it. Such further complication in extracting the gluon shadowing from 
the BRAHMS data is, however, not considered here, since the inclusion of the charge 
separation becomes more reliable only in NLO. 

Our next goal is to perform this analysis in NLO, as well as, when the data be- 
come finalized, include other RHIC data sets, such as photon production in d-l-Au and 
Au-|-Au from PHENIX, into the analysis. In general, any further constraints for the 
gluon distributions are more than welcome. For example, the ratio -RdAu for D mesons 
to be (hopefully soon) measured at RHIC will be extremely useful, at any rapidity. The 
cleanest environment for the nPDFs measurements would be in the DIS experiments 
at eRHIC [16] and LHeC colliders now being discussed. Before the possible realization 
of these machines, we hope that the present study in its part demonstrates how impor- 
tant it would be for the correct determination of universal nPDFs to have a systematic 
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proton-nucleus program also at the LHC: further constraints for nuclear gluons in the 
yet unexplored regions of the x, plane are absolutely necessary for understanding 
QCD parton dynamics in high-energy nuclear and hadronic collisions. 
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